Below you can find a link to the MATLAB code calculating the Galerkin advection-diffusion problem for the Eriksson-Johnson model problem.

Let us recall first the advection-diffusion problem, the Eriksson-Johnson model problem (described, for example, in [1]) defined on a square area \( \Omega = [0,1]^2 \) in the following way: We see \( u \) such as
\( a(u,v)=l(v) \forall v \) where
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) is the extension of the Dirichlet boundary condition into the entire area, and
\( \beta = (1,0) \) represents the wind velocity blowing from the left side, and \( \epsilon = 10^{-2} \) denotes the diffusion coefficient. The quantity \( Pe=1/ \epsilon = 100 \) is called the Peclet number, and it defines the numerical sensitivity of the advection-diffusion problem. The definition of the Peclet number presented here is valid for unit wind speed \( \| \beta \|_2=\sqrt{\beta_x^2+\beta^2_y}=1 \). In the general case, the Peclet number defined on the individual grid elements according to the formula is considered
\( Pe=\| \beta \|_2*h/(2*\epsilon) \). The Peclet number defined in this way can be used to test the stability of the simulation, the Peclet number less than one gives us the condition for the stability of the simulation on the element.

In line 835 we create the knot vector with 5 elements for quadratic polynomials
\( knot = simple\_knot(5, 2); \)
In line 836, the diffusion coefficient is set (for an advection coefficient of 1).
\( epsilon = 2e-2; \)

The code execution is also possible in the free Octave environment.
Download code or see ( Appendix 3A ).

The code is run by opening it in Octave and typing a command
\( advection' \).
After a moment of calculation, the code opens an additional window and draws a numerical and exact solution to it.
The code also computes the norm L2 and H1 from the solution
\( Error: L2 24.00 procent, H1 72.08 procent \)
and compares it to the L2 and H1 standards from the projection of the exact solution onto the space of B-spline functions spanning a given grid. \( Best possible: L2 1.56 procent, H1 53.16 procent \)
Note that this norm will not be equal to zero, because the exact solution is not given by the B-spline functions, and even if we know the exact solution, it is not given by the B-spline functions, and its projection on the B-spline function space is also an approximation of the exact solution.